{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 60,
   "id": "8017a454-1988-49d9-b964-ae31ebe24ced",
   "metadata": {},
   "outputs": [],
   "source": [
    "import struct\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import xarray as xr\n",
    "import matplotlib.pyplot as plt\n",
    "import fiona\n",
    "from matplotlib.path import Path\n",
    "import matplotlib as mpl"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "a3c68087-0f07-4555-9a13-712fba9bfed6",
   "metadata": {},
   "outputs": [],
   "source": [
    "###function to calculate distance between points\n",
    "\n",
    "def great_circle_dist_moon(p1,p2,r=1737.0):\n",
    "    \"\"\"\n",
    "    Calculate great-circle distance between two points\n",
    "    \n",
    "    Input:\n",
    "        p1 = [lon, lat] \n",
    "        p2 = [lon, lat]\n",
    "\n",
    "        r = Radius (Default 1737.0 km) \n",
    "    Output:\n",
    "        d = distance (Default km)\n",
    "    \"\"\"\n",
    "    \n",
    "    lon1 = np.deg2rad(p1[0])\n",
    "    lon2 = np.deg2rad(p2[0])\n",
    "    lat1 = np.deg2rad(p1[1])\n",
    "    lat2 = np.deg2rad(p2[1])\n",
    "    \n",
    "    cent_ang = np.arccos(np.sin(lat1)*np.sin(lat2) + np.cos(lat1)*np.cos(lat2)*np.cos(np.abs(lon1-lon2)))\n",
    "    d = r*cent_ang\n",
    "    return d"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "710b5a6d-cd76-4198-8184-cbe131069033",
   "metadata": {
    "tags": []
   },
   "source": [
    "# Load Datasets:\n",
    "\n",
    "### Thorium Data\n",
    "https://pds-geosciences.wustl.edu/missions/lunarp/reduced/thoriumhd.txt "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "b8ab003a-48a7-4b0e-b2a4-cf707d4a4b6f",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/var/folders/zj/46zkt8ns517_676qkc7f89tm0000gv/T/ipykernel_46798/3203252587.py:2: ParserWarning: Falling back to the 'python' engine because the 'c' engine does not support regex separators (separators > 1 char and different from '\\s+' are interpreted as regex); you can avoid this warning by specifying engine='python'.\n",
      "  dat_th = pd.read_csv('../data/thoriumhd.txt', skiprows = 13, sep=',\\s+', lineterminator='\\n', names=[\"lat_min\",\"lat_max\",\"lon_min\",\"lon_max\",\"ppm\"])\n"
     ]
    }
   ],
   "source": [
    "## Load thorium data\n",
    "dat_th = pd.read_csv('../data/thoriumhd.txt', skiprows = 13, sep=',\\s+', lineterminator='\\n', names=[\"lat_min\",\"lat_max\",\"lon_min\",\"lon_max\",\"ppm\"])\n",
    "\n",
    "ppm = dat_th.ppm.values.reshape([360,720])\n",
    "thlon = dat_th.lon_min.values.reshape([360,720])\n",
    "thlat = dat_th.lat_min.values.reshape([360,720])\n",
    "\n",
    "\n",
    "## transform into xarray Dataset\n",
    "thor = xr.Dataset(dict(th_ppm = ([\"x\",\"y\"],ppm)),coords = dict(lat=([\"x\",\"y\"],thlat),lon=([\"x\",\"y\"],thlon))) \n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "dc73fbb4-e617-4095-80d2-87fa684537c6",
   "metadata": {},
   "source": [
    "### Crater Survey \n",
    "\n",
    "https://astrogeology.usgs.gov/search/map/Moon/Research/Craters/lunar_crater_database_robbins_2018 "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "b652ace9-9682-45aa-9802-c8af71dbbdc6",
   "metadata": {},
   "outputs": [],
   "source": [
    "#Load Robbins 2018 crater survey  \n",
    "df = pd.read_csv('../data/lunar_crater_database_robbins_2018.csv')\n",
    "df = df.set_index('CRATER_ID')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b2721fa3-d7c4-4362-ac21-86f18bb0befb",
   "metadata": {},
   "source": [
    "## Generate Mare Mask\n",
    "\n",
    "https://wms.lroc.asu.edu/lroc/view_rdr/SHAPEFILE_LROC_GLOBAL_MARE"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f26b5a7f-a7ab-4b8e-b05b-dc6a83c259fd",
   "metadata": {
    "tags": []
   },
   "outputs": [],
   "source": [
    "shape180 = fiona.open('../data/LROC_GLOBAL_MARE_180.SHP') #open shapefile\n",
    "\n",
    "\n",
    "plon, plat = thlon.flatten(), thlat.flatten()\n",
    "points = np.vstack((plon,plat)).T \n",
    "th_mask_180 = np.zeros((360,720), dtype = bool)\n",
    "for i in range(len(shape180)):\n",
    "    try:\n",
    "        perim = np.array((shape180[i]['geometry']['coordinates'][0])) #unpack coordinates\n",
    "        p = Path(perim) # make a polygon\n",
    "        grid = p.contains_points(points)\n",
    "        local_mask = grid.reshape(360,720)\n",
    "        th_mask_180 = th_mask_180 + local_mask\n",
    "\n",
    "\n",
    "    #some are multi polygons for some reason. anyway\n",
    "    except ValueError:\n",
    "        for j in range(len(shape180[i]['geometry']['coordinates'])):\n",
    "            perim = np.array((shape180[i]['geometry']['coordinates'][j][0])) \n",
    "        \n",
    "            p = Path(perim) # make a polygon\n",
    "            grid = p.contains_points(points)\n",
    "            local_mask = grid.reshape(360,720)\n",
    "            th_mask_180 = th_mask_180 + local_mask\n",
    "np.savetxt(\"mare_mask_0.5deg.txt\", th_mask_180, fmt=\"%5i\")\n",
    "         "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d8438970-5886-4cc1-9be6-8e1ba96ce660",
   "metadata": {},
   "source": [
    "### or load Mare Mask"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e245cf49-25b0-49ed-bd77-0a2d4aa7be4f",
   "metadata": {},
   "outputs": [],
   "source": [
    "th_mask_180 = np.loadtxt(\"mare_mask_0.5deg.txt\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "18b3f886-a105-438b-a9af-dd3ed541c339",
   "metadata": {},
   "source": [
    "### Apply Mare Mask to Thorium Data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9e0a0905-898b-497f-b68d-608cb34d05f3",
   "metadata": {},
   "outputs": [],
   "source": [
    "th_mask_180 = np.loadtxt(\"mare_mask_0.5deg.txt\")\n",
    "thor_nomare = thor.where(th_mask_180 == False)\n",
    "weight = np.cos(np.deg2rad(np.abs(thor.lat)))\n",
    "\n",
    "\n",
    "thor_nomare['weights'] = (('x','y'),weight.values)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f1017c5a-b2fe-45a7-939e-427f2c869cb8",
   "metadata": {
    "tags": []
   },
   "source": [
    "# Constants and Functions\n",
    "\n",
    "### LP GRS Instrument Response Function \n",
    "(Lawrence et al. 2003, Kivelson 1995)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "id": "0188174a-8abd-4fa1-b936-18d482aa3abd",
   "metadata": {},
   "outputs": [],
   "source": [
    "A = 1; \n",
    "\n",
    "sigma = lambda h: 0.704*h + 1.39\n",
    "kappa = lambda h: -4.87*(10**(-4))*h + 0.631\n",
    "\n",
    "w = lambda D, h: A*(1 + D**2/(2*sigma(h)**2))**(-kappa(h)-1)\n",
    "\n",
    "D = np.linspace(-600,600,100) #Distance from sub-nadir point (km)\n",
    "h = 30 #Altitude (km) \n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b532ab01-03eb-404d-8935-7b8ee9a48daa",
   "metadata": {},
   "source": [
    "# Calculate TSTAs\n",
    "\n",
    "### Select crater subset by diameter:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "id": "4f6abf8a-ec21-4c3b-a75d-af7713b14108",
   "metadata": {},
   "outputs": [],
   "source": [
    "### select diameters between 20 and 200 km\n",
    "dfsubset = df.loc[(df['DIAM_CIRC_IMG'] >= 20.0)&(df['DIAM_CIRC_IMG'] <= 200.0)]"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ded54dfc-7d41-4f06-a703-61bea0a6069c",
   "metadata": {},
   "source": [
    "### Compute: \n",
    "\n",
    "Expected run-time of 1 - 2 hours"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "id": "b50864d8-3e7f-4f64-b7e8-acce369034d1",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<==|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||•==> 99%\r"
     ]
    }
   ],
   "source": [
    "th_20_200_circ = []\n",
    "tot = dfsubset.shape[0]\n",
    "\n",
    "for p in range(tot):\n",
    "    basin = dfsubset.iloc[p]\n",
    "    lat = basin.LAT_CIRC_IMG\n",
    "    lon = basin.LON_CIRC_IMG\n",
    "    dc = basin.DIAM_CIRC_IMG\n",
    "        \n",
    "    ## correcting to -180 to 180 longitudes\n",
    "    if lon > 180.0:\n",
    "        lon_th = lon - 360.0\n",
    "    else:\n",
    "        lon_th = lon\n",
    "        \n",
    "    rc = dc/2 #Radius of crater\n",
    "\n",
    "    r_ej = 2*rc #Radius of continuos ejecta (Melosh 1989)\n",
    "    \n",
    "    r_back = r_ej + 300 ###Radius of background region\n",
    "    \n",
    "    \n",
    "    ### Generate masks for each region and data averages\n",
    "    \n",
    "    ### distances from center of crater\n",
    "    gc_dist = great_circle_dist_moon([thor.lon.values,thor.lat.values], [lon_th,lat])\n",
    "    \n",
    "    ### mask of crater interior\n",
    "    cratermaskarr = np.zeros_like(gc_dist)\n",
    "    cratermaskarr[np.where(gc_dist < rc)] = 1 \n",
    "\n",
    "\n",
    "    ### mask of ejecta annulus\n",
    "    ejectamaskarr = np.zeros_like(gc_dist)\n",
    "    ejectamaskarr[np.where((gc_dist < r_ej ) & (gc_dist > rc))] = 1 \n",
    "    \n",
    "    ### mask of background annulus\n",
    "    backmaskarr = np.zeros_like(gc_dist)\n",
    "    backmaskarr[np.where((gc_dist < r_back) & (gc_dist > r_ej))] = 1 \n",
    "\n",
    "    #### Compute area-weighted mean (excluding mare pixels)\n",
    "\n",
    "    crat_th_m = thor_nomare.th_ppm.where(cratermaskarr == 1).weighted(weight).mean()\n",
    "    \n",
    "    ej_th_m = thor_nomare.th_ppm.where(ejectamaskarr == 1).weighted(weight).mean()\n",
    "    \n",
    "    back_m = thor_nomare.th_ppm.where(backmaskarr == 1).weighted(weight).mean()\n",
    "\n",
    "\n",
    "    ##For each region, if there are more than 2 non-mare pixels, compute area-weighted standard deviation\n",
    "    \n",
    "    M_crater = np.count_nonzero(cratermaskarr)\n",
    "\n",
    "    if M_crater > 1:\n",
    "        std_w_crater = np.sqrt((weight*(thor_nomare.th_ppm.where(cratermaskarr == 1) - crat_th_m.item())**2).sum()/\n",
    "                        (((M_crater-1)/M_crater)*weight.where(cratermaskarr==1).sum())).item()\n",
    "    else:\n",
    "        std_w_crater = np.nan\n",
    "\n",
    "    M_ej = np.count_nonzero(ejectamaskarr)\n",
    "    if M_ej > 1:\n",
    "        std_w_ej = np.sqrt((weight*(thor_nomare.th_ppm.where(ejectamaskarr == 1) - ej_th_m.item())**2).sum()/\n",
    "                        (((M_ej-1)/M_ej)*weight.where(ejectamaskarr==1).sum())).item()\n",
    "    else:\n",
    "        std_w_ej = np.nan\n",
    "    M_back = np.count_nonzero(backmaskarr)\n",
    "\n",
    "    if M_back > 1:\n",
    "        std_w_back = np.sqrt((weight*(thor_nomare.th_ppm.where(backmaskarr == 1) - back_m.item())**2).sum()/\n",
    "                        (((M_back-1)/M_back)*weight.where(backmaskarr==1).sum())).item()\n",
    "    else:\n",
    "        std_w_back = np.nan\n",
    "\n",
    "\n",
    "    ###### Radial averaging \n",
    "    \n",
    "    ###mask out everything outside of the three regions (everything greater than the background region) \n",
    "    \n",
    "    relevant = np.zeros_like(gc_dist)\n",
    "    relevant[np.where(gc_dist< r_back)] = 1.0\n",
    "    distarr = gc_dist * relevant\n",
    "    \n",
    "    thor_nomare['distarr'] = (('x','y'),distarr) #distance up to background border\n",
    "    \n",
    "    \n",
    "    bins = np.linspace(0,r_back,100) #Bins for radial averaging\n",
    "    \n",
    "    gb_all = thor_nomare.groupby_bins(thor_nomare.distarr,bins[1:]) #separate into bins by distance\n",
    "    \n",
    "    ngroups = len(gb_all.groups)\n",
    "\n",
    "    prof = np.zeros((ngroups,2))\n",
    "\n",
    "    i = 0\n",
    "    #weighted average of each bin\n",
    "    for gr in gb_all.groups:\n",
    "        keep = ~np.isnan(gb_all[gr].th_ppm.values) # non nan values\n",
    "        if keep.any() == True:\n",
    "            prof[i,:] = [gr.left,np.average(gb_all[gr].th_ppm.values[keep],weights = gb_all[gr].weights.values[keep])]\n",
    "        else:\n",
    "            prof[i,:] = [gr.left,np.nan]\n",
    "        #print(gr, np.average(gb_all[gr].th_ppm.values,weights = gb_all[gr].weights.values))\n",
    "        i = i + 1\n",
    "    sort_prof = prof[prof[:, 0].argsort()] #sort back into order by distance from center\n",
    "\n",
    "    #set up for convolution by reflecting over y axis \n",
    "    doubled = np.append(np.flip(sort_prof[:,1]),sort_prof[:,1])\n",
    "    D_prof = np.append((-1)*np.flip(sort_prof[:,0]),sort_prof[:,0])\n",
    "\n",
    "    #mask out distances > distance array\n",
    "    \n",
    "    double_mask = doubled.copy()\n",
    "    double_mask[np.where(np.abs(D_prof) > r_back)] = 0 \n",
    "\n",
    "    #### convolve with Instrument Response Function (IRF)\n",
    "    convdat = np.convolve(double_mask, w(D,h))\n",
    "    D_conv = np.linspace(-600,600,len(convdat))\n",
    "    \n",
    "    ### Compute Mean Observed +Convolved values\n",
    "\n",
    "\n",
    "    O3 = np.zeros((3,1))\n",
    "    O3[0] = np.mean(convdat[np.where(np.abs(D_conv)  < rc)])\n",
    "    O3[1] = np.mean(convdat[np.where((np.abs(D_conv)  < r_ej) & (np.abs(D_conv) > rc))])\n",
    "    O3[2] = np.mean(convdat[np.where((np.abs(D_conv) > r_ej)&(np.abs(D_conv) < r_back))])\n",
    "    \n",
    "    \n",
    "    \n",
    "\n",
    "\n",
    "    ### Make masks where each region is unity and everything else is 0 and convolve with IRF\n",
    "    \n",
    "    #Crater Interior\n",
    "    modc = np.zeros_like(D_prof)\n",
    "    modc[np.where(np.abs(D_prof) < rc)] = 1\n",
    "\n",
    "    convc = np.convolve(modc,w(D,h))\n",
    "\n",
    "\n",
    "\n",
    "    #Crater Ejecta\n",
    "    mode = np.zeros_like(D_prof)\n",
    "    mode[np.where((np.abs(D_prof) > rc)&(np.abs(D_prof) < r_ej))] = 1\n",
    "\n",
    "    conve = np.convolve(mode,w(D,h))\n",
    "\n",
    "\n",
    "    #Background region\n",
    "    modx = np.zeros_like(D_prof)\n",
    "    modx[np.where((np.abs(D_prof) > r_ej)&(np.abs(D_prof) < r_back))] = 1\n",
    "\n",
    "    convx = np.convolve(modx,w(D,h))\n",
    "    \n",
    "    \n",
    "\n",
    "    ## Set up Instrument Response Matrix for deconvolution (how much would a signal of 1 from each region contribute to the signal observed in each region)\n",
    "\n",
    "    C3 = np.zeros((3,3))\n",
    "    C3[0,0] = np.mean(convc[np.where(np.abs(D_conv)  < rc)]) #cii\n",
    "    C3[0,1] = np.mean(conve[np.where(np.abs(D_conv)  < rc)]) #cie\n",
    "    C3[0,2] = np.mean(convx[np.where(np.abs(D_conv)  < rc)]) #cix\n",
    "\n",
    "    C3[1,0] = np.mean(convc[np.where((np.abs(D_conv)  < r_ej) & (np.abs(D_conv) > rc))]) #cei\n",
    "    C3[1,1] = np.mean(conve[np.where((np.abs(D_conv)  < r_ej) & (np.abs(D_conv) > rc))]) #cee\n",
    "    C3[1,2] = np.mean(convx[np.where((np.abs(D_conv)  < r_ej) & (np.abs(D_conv) > rc))]) #cex\n",
    "\n",
    "    C3[2,0] = np.mean(convc[np.where((np.abs(D_conv) > r_ej)&(np.abs(D_conv) < r_back))]) #cxi\n",
    "    C3[2,1] = np.mean(conve[np.where((np.abs(D_conv) > r_ej)&(np.abs(D_conv) < r_back))]) #cxe\n",
    "    C3[2,2] = np.mean(convx[np.where((np.abs(D_conv) > r_ej)&(np.abs(D_conv) < r_back))]) #cxx\n",
    "    \n",
    "    #deconvolve to find \"True Surface Thorium Anomalies\" (TSTA)\n",
    "    try:\n",
    "        T3 = np.matmul(np.linalg.inv(C3),O3)\n",
    "        boundarr = gc_dist - rc #\n",
    "\n",
    "        out_idx = (boundarr < 5) & (boundarr > -5)\n",
    "        #build dictionary for pandas array\n",
    "\n",
    "        d = {\n",
    "            'idx' : basin.name,\n",
    "            'lat' : lat,\n",
    "            'lon' : lon_th,\n",
    "            'rc' : rc,\n",
    "            'crater_m' : crat_th_m.values.item(),\n",
    "            'crater_std': std_w_crater,\n",
    "            'ejecta_m': ej_th_m.values.item(),\n",
    "            'ejecta_std' : std_w_ej,\n",
    "            'back_m' : back_m.values.item(),\n",
    "            'back_std': std_w_back,\n",
    "            'floor_m_t' : T3[0].item(),\n",
    "            'ejecta_m_t': T3[1].item(),\n",
    "            'back_m_t': T3[2].item(),\n",
    "            'lon_perimeter' : thor.lon.values[out_idx],  # coordinates of crater perimeter\n",
    "            'lat_perimeter' : thor.lat.values[out_idx]\n",
    "        }\n",
    "        th_20_200_circ.append(d)\n",
    "    except np.linalg.LinAlgError as err:\n",
    "        print('Singular matrix', end=\"\\r\")\n",
    "\n",
    "    ###loading bar \n",
    "    if (p > 0) & (p % (tot//100) == 0): \n",
    "        done = int(p/(tot//100 )) -1\n",
    "        togo = int(100-done)\n",
    "        print(\"<==\" + done*'|'+ togo*'•' + \"==> \" + str(done) + \"%\", end=\"\\r\")\n",
    "print(\"<==\" + 100*'|'+ \"==> \" + str(done) + \"%\", end=\"\\r\")\n",
    "##make pandas dataframe\n",
    "th_abun_circ = pd.DataFrame(th_20_200_circ)\n",
    "std_pd_circ = th_abun_circ.set_index('idx')\n",
    "### save as csv\n",
    "std_pd_circ.to_csv('../results/tsta_20_200.csv')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "5a255521-189d-4909-bfb0-84edbf0041c4",
   "metadata": {},
   "source": [
    "# Caclulate TSTAs for Randomized Crater Locations\n",
    "\n",
    "Same procedure but with randomized center coordinates\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a983ec0d-42a1-4673-9530-dc324c965de9",
   "metadata": {},
   "outputs": [],
   "source": [
    "import random \n",
    "import warnings\n",
    "\n",
    "warnings.filterwarnings(\"ignore\", message=\"Degrees of freedom\")\n",
    "\n",
    "###"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "id": "a544a6ec-0e6b-463d-9fb2-95799a78e894",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "<==|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||•==> 99%\r"
     ]
    }
   ],
   "source": [
    "import random\n",
    "th_random_set = []\n",
    "tot = dfsubset.shape[0]\n",
    "\n",
    "lats = np.linspace(-90,90,720)\n",
    "lons = np.linspace(0,360,1440) \n",
    "for p in range(tot):\n",
    "    basin = dfsubset.iloc[p]\n",
    "    lat = np.random.choice(lats)\n",
    "    lon = np.random.choice(lons)\n",
    "    dc = basin.DIAM_CIRC_IMG\n",
    "        \n",
    "    ## correcting to -180 to 180 longitudes\n",
    "    if lon > 180.0:\n",
    "        lon_th = lon - 360.0\n",
    "    else:\n",
    "        lon_th = lon\n",
    "        \n",
    "    rc = dc/2 #Radius of crater\n",
    "\n",
    "    r_ej = 2*rc #Radius of continuos ejecta (Melosh 1989)\n",
    "    \n",
    "    r_back = r_ej + 300 ###Radius of background region\n",
    "    \n",
    "    \n",
    "    ### Generate masks for each region and data averages\n",
    "    \n",
    "    ### distances from center of crater\n",
    "    gc_dist = great_circle_dist_moon([thor.lon.values,thor.lat.values], [lon_th,lat])\n",
    "    \n",
    "    ### mask of crater interior\n",
    "    cratermaskarr = np.zeros_like(gc_dist)\n",
    "    cratermaskarr[np.where(gc_dist < rc)] = 1 \n",
    "\n",
    "\n",
    "    ### mask of ejecta annulus\n",
    "    ejectamaskarr = np.zeros_like(gc_dist)\n",
    "    ejectamaskarr[np.where((gc_dist < r_ej ) & (gc_dist > rc))] = 1 \n",
    "    \n",
    "    ### mask of background annulus\n",
    "    backmaskarr = np.zeros_like(gc_dist)\n",
    "    backmaskarr[np.where((gc_dist < r_back) & (gc_dist > r_ej))] = 1 \n",
    "\n",
    "    #### Compute area-weighted mean (excluding mare pixels)\n",
    "\n",
    "    crat_th_m = thor_nomare.th_ppm.where(cratermaskarr == 1).weighted(weight).mean()\n",
    "    \n",
    "    ej_th_m = thor_nomare.th_ppm.where(ejectamaskarr == 1).weighted(weight).mean()\n",
    "    \n",
    "    back_m = thor_nomare.th_ppm.where(backmaskarr == 1).weighted(weight).mean()\n",
    "\n",
    "\n",
    "    ##For each region, if there are more than 2 non-mare pixels, compute area-weighted standard deviation\n",
    "    \n",
    "    M_crater = np.count_nonzero(cratermaskarr)\n",
    "\n",
    "    if M_crater > 1:\n",
    "        std_w_crater = np.sqrt((weight*(thor_nomare.th_ppm.where(cratermaskarr == 1) - crat_th_m.item())**2).sum()/\n",
    "                        (((M_crater-1)/M_crater)*weight.where(cratermaskarr==1).sum())).item()\n",
    "    else:\n",
    "        std_w_crater = np.nan\n",
    "\n",
    "    M_ej = np.count_nonzero(ejectamaskarr)\n",
    "    if M_ej > 1:\n",
    "        std_w_ej = np.sqrt((weight*(thor_nomare.th_ppm.where(ejectamaskarr == 1) - ej_th_m.item())**2).sum()/\n",
    "                        (((M_ej-1)/M_ej)*weight.where(ejectamaskarr==1).sum())).item()\n",
    "    else:\n",
    "        std_w_ej = np.nan\n",
    "    M_back = np.count_nonzero(backmaskarr)\n",
    "\n",
    "    if M_back > 1:\n",
    "        std_w_back = np.sqrt((weight*(thor_nomare.th_ppm.where(backmaskarr == 1) - back_m.item())**2).sum()/\n",
    "                        (((M_back-1)/M_back)*weight.where(backmaskarr==1).sum())).item()\n",
    "    else:\n",
    "        std_w_back = np.nan\n",
    "\n",
    "\n",
    "    ###### Radial averaging \n",
    "    \n",
    "    ###mask out everything outside of the three regions (everything greater than the background region) \n",
    "    \n",
    "    relevant = np.zeros_like(gc_dist)\n",
    "    relevant[np.where(gc_dist< r_back)] = 1.0\n",
    "    distarr = gc_dist * relevant\n",
    "    \n",
    "    thor_nomare['distarr'] = (('x','y'),distarr) #distance up to background border\n",
    "    \n",
    "    \n",
    "    bins = np.linspace(0,r_back,100) #Bins for radial averaging\n",
    "    \n",
    "    gb_all = thor_nomare.groupby_bins(thor_nomare.distarr,bins[1:]) #separate into bins by distance\n",
    "    \n",
    "    ngroups = len(gb_all.groups)\n",
    "\n",
    "    prof = np.zeros((ngroups,2))\n",
    "\n",
    "    i = 0\n",
    "    #weighted average of each bin\n",
    "    for gr in gb_all.groups:\n",
    "        keep = ~np.isnan(gb_all[gr].th_ppm.values) # non nan values\n",
    "        if keep.any() == True:\n",
    "            prof[i,:] = [gr.left,np.average(gb_all[gr].th_ppm.values[keep],weights = gb_all[gr].weights.values[keep])]\n",
    "        else:\n",
    "            prof[i,:] = [gr.left,np.nan]\n",
    "        #print(gr, np.average(gb_all[gr].th_ppm.values,weights = gb_all[gr].weights.values))\n",
    "        i = i + 1\n",
    "    sort_prof = prof[prof[:, 0].argsort()] #sort back into order by distance from center\n",
    "\n",
    "    #set up for convolution by reflecting over y axis \n",
    "    doubled = np.append(np.flip(sort_prof[:,1]),sort_prof[:,1])\n",
    "    D_prof = np.append((-1)*np.flip(sort_prof[:,0]),sort_prof[:,0])\n",
    "\n",
    "    #mask out distances > distance array\n",
    "    \n",
    "    double_mask = doubled.copy()\n",
    "    double_mask[np.where(np.abs(D_prof) > r_back)] = 0 \n",
    "\n",
    "    #### convolve with Instrument Response Function (IRF)\n",
    "    convdat = np.convolve(double_mask, w(D,h))\n",
    "    D_conv = np.linspace(-600,600,len(convdat))\n",
    "    \n",
    "    ### Compute Mean Observed +Convolved values\n",
    "\n",
    "\n",
    "    O3 = np.zeros((3,1))\n",
    "    O3[0] = np.mean(convdat[np.where(np.abs(D_conv)  < rc)])\n",
    "    O3[1] = np.mean(convdat[np.where((np.abs(D_conv)  < r_ej) & (np.abs(D_conv) > rc))])\n",
    "    O3[2] = np.mean(convdat[np.where((np.abs(D_conv) > r_ej)&(np.abs(D_conv) < r_back))])\n",
    "    \n",
    "    \n",
    "    \n",
    "\n",
    "\n",
    "    ### Make masks where each region is unity and everything else is 0 and convolve with IRF\n",
    "    \n",
    "    #Crater Interior\n",
    "    modc = np.zeros_like(D_prof)\n",
    "    modc[np.where(np.abs(D_prof) < rc)] = 1\n",
    "\n",
    "    convc = np.convolve(modc,w(D,h))\n",
    "\n",
    "\n",
    "\n",
    "    #Crater Ejecta\n",
    "    mode = np.zeros_like(D_prof)\n",
    "    mode[np.where((np.abs(D_prof) > rc)&(np.abs(D_prof) < r_ej))] = 1\n",
    "\n",
    "    conve = np.convolve(mode,w(D,h))\n",
    "\n",
    "\n",
    "    #Background region\n",
    "    modx = np.zeros_like(D_prof)\n",
    "    modx[np.where((np.abs(D_prof) > r_ej)&(np.abs(D_prof) < r_back))] = 1\n",
    "\n",
    "    convx = np.convolve(modx,w(D,h))\n",
    "    \n",
    "    \n",
    "\n",
    "    ## Set up Instrument Response Matrix for deconvolution (how much would a signal of 1 from each region contribute to the signal observed in each region)\n",
    "\n",
    "    C3 = np.zeros((3,3))\n",
    "    C3[0,0] = np.mean(convc[np.where(np.abs(D_conv)  < rc)]) #cii\n",
    "    C3[0,1] = np.mean(conve[np.where(np.abs(D_conv)  < rc)]) #cie\n",
    "    C3[0,2] = np.mean(convx[np.where(np.abs(D_conv)  < rc)]) #cix\n",
    "\n",
    "    C3[1,0] = np.mean(convc[np.where((np.abs(D_conv)  < r_ej) & (np.abs(D_conv) > rc))]) #cei\n",
    "    C3[1,1] = np.mean(conve[np.where((np.abs(D_conv)  < r_ej) & (np.abs(D_conv) > rc))]) #cee\n",
    "    C3[1,2] = np.mean(convx[np.where((np.abs(D_conv)  < r_ej) & (np.abs(D_conv) > rc))]) #cex\n",
    "\n",
    "    C3[2,0] = np.mean(convc[np.where((np.abs(D_conv) > r_ej)&(np.abs(D_conv) < r_back))]) #cxi\n",
    "    C3[2,1] = np.mean(conve[np.where((np.abs(D_conv) > r_ej)&(np.abs(D_conv) < r_back))]) #cxe\n",
    "    C3[2,2] = np.mean(convx[np.where((np.abs(D_conv) > r_ej)&(np.abs(D_conv) < r_back))]) #cxx\n",
    "    \n",
    "    #deconvolve to find \"True Surface Thorium Anomalies\" (TSTA)\n",
    "    try:\n",
    "        T3 = np.matmul(np.linalg.inv(C3),O3)\n",
    "        boundarr = gc_dist - rc #\n",
    "\n",
    "        out_idx = (boundarr < 5) & (boundarr > -5)\n",
    "        #build dictionary for pandas array\n",
    "\n",
    "        d = {\n",
    "            'idx' : basin.name,\n",
    "            'lat' : lat,\n",
    "            'lon' : lon_th,\n",
    "            'rc' : rc,\n",
    "            'crater_m' : crat_th_m.values.item(),\n",
    "            'crater_std': std_w_crater,\n",
    "            'ejecta_m': ej_th_m.values.item(),\n",
    "            'ejecta_std' : std_w_ej,\n",
    "            'back_m' : back_m.values.item(),\n",
    "            'back_std': std_w_back,\n",
    "            'floor_m_t' : T3[0].item(),\n",
    "            'ejecta_m_t': T3[1].item(),\n",
    "            'back_m_t': T3[2].item(),\n",
    "            'lon_perimeter' : thor.lon.values[out_idx],  # coordinates of crater perimeter\n",
    "            'lat_perimeter' : thor.lat.values[out_idx]\n",
    "        }\n",
    "        th_random.append(d)\n",
    "    except np.linalg.LinAlgError as err:\n",
    "        print('Singular matrix', end=\"\\r\")\n",
    "\n",
    "    ###loading bar \n",
    "    if (p > 0) & (p % (tot//100) == 0): \n",
    "        done = int(p/(tot//100 )) -1\n",
    "        togo = int(100-done)\n",
    "        print(\"<==\" + done*'|'+ togo*'•' + \"==> \" + str(done) + \"%\", end=\"\\r\")\n",
    "print(\"<==\" + 100*'|'+ \"==> \" + str(done) + \"%\", end=\"\\r\")\n",
    "##make pandas dataframe\n",
    "th_abun_random = pd.DataFrame(th_random_set)\n",
    "thabun = th_abun_random.set_index('idx')\n",
    "### save as csv\n",
    "thabun.to_csv('../results/tsta_random.csv')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6554b09f-ac29-4c52-ae64-2b8132bed019",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
